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ABSTRACT 


Failure mechanisms of materials under very high strains experienced 
at and ahead of the crack tip such as the formation, growth and interac- 
tion of microvoids in ductile materials, microcracks in brittle solids or 
crazes in polymers and adhesives are represented by one-dimensional , non- 
linear stress-strain relations possessing different post-yield softening 
(strain-softening) behavior. These reflect different ways by which the 
material loses capacity to carry load up to fracture or total seperation. 
A DCB type specimen is considered in this study. 1 2 The nonlinear material 
is confined to a thin strip between the two elastic beams loaded by a 
wedge. The problem is first modelled as a beam on a nonlinear foundation. 
The pertinent equation is solved numerically as a two-point boundary value 
problem for both the stationary and the quasi-statically propagating 
crack. A finite element model is then used to model the problem in more 
detail in order to assess the adequacy of the beam model for the reduction 
of experimental data to determine in-situ properties of the thin inter- 
layer. 

It is found that the energy release rate derived by assuming built-in 
conditions at the crack tip can be used to calculate the fracture (sur- 
face) energy more accurately and conveniently than Berry's scheme [2] even 
in cases where the built-in assumption is apparently invalid. The ana- 
lyses suggest that it is possible to infer from detailed macroscopic meas- 
urements of the deformations of the beam prior to and during crack growth 
the approximate characteristics of the complete (uniaxial) material 
stress-strain behavior of the cohesive interlayer including loading and 
strain-sof tening . ® 


1. This paper is a shortened version of reference 1. 

2. An experimental program is being conducted in our laboratories to study 
the stress-strain characteristics of a number of polymeric solids using 
the DCB model discussed in this paper. 
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1. INTRODUCTION 

The double cantilever beam (DCB) specimen has been used extensively 
in crack propagation studies due to its simple geometry which is attrac- 
tive from both experimental and theoretical standpoints. Berry [3] inves- 
tigated the implications of the Griffith fracture criterion for the DCB 
specimen and showed that for the built-in beam, fracture initiated when 
the moment at the crack tip reached a critical value. Bilek and Burns [4] 
in their dynamic analysis derived the same equations through the applica- 
tion of Hamilton's principle. Similar results were also obtained by 
Steverding and Lehnigk [5] using a different approach. Kanninen [6] 
employed the model of a beam on an elastic foundation to relax the built- 
in constraint and later [7] extended his analysis to the dynamic case 
introducing a Timoshenko beam on a generalized elastic foundation model to 
account for shear deformation and rotary inertia. 

In the present work, we generalize the foundation to a nonlinear one 
characterized by an initially linear elastic stress-strain relation and an 
unloading tail reflecting loss of load carrying ability (strain- 
softening). Subsequently, fracture occurs at some critical strain at 
which the foundation stiffness drops to zero. Thus the model allows the 
crack to propagate without any additional prescription of a failure cri- 
terion, so that the nonlinear DCB analysis becomes a potential tool for 
determining the in-situ, nonlinear material characteristics of the bond- 
interlayer. As complete stress-strain characterizations of real materials 
are not yet available, we resort first to some idealized (hypothetical) 
material models in an attempt to extract from the solution of the problem 
certain measurable macroscopic quantities that might allow one to charac- 
terize the complete material behavior in the continuum sense: In this we 
assume that a sufficiently large number of microvoids, microcracks or 
craze fibrils are present in each small volume element under load at the 
"crack front" such that the damage-induced loss of load carrying capacity 
can be meaningfully averaged and represented as a continuum response. 
Studies involving porous and damaged materials are numerous, for instance 
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McClintock [8] , Gurson [9] who developed approximate yield criteria and 
flow rules for dilatant ductile materials; Needleman [10] , Berg [11] who 
proposed a continuum model for plastic deformation of microporous metal 
aggregates; Rice and Needleman [12] studied metals with constitutive 
dependence on hydrostatic stress (which promotes void nucleation and 
growth) in sheet metal forming processes; Krajcinovic [13] presented a 
constitutive model for material containing flat planar microcracks such as 
concrete; see also Dougill [14], Rudnicki [15] and Bazant [16]. In our 
studies we are primarily interested in the nonlinear response of homogene- 
ous or multiphase polymers such as those used in advanced composites and 
adhesives. 

To simulate crack propagation in our finite element study, we employ 
a scheme different from the ones proposed thus far. The usual node 
release methods proposed by Andersson [17], Kfouri & Miller [18], Rydholm 
et al . [19] and Malluck & King [20, 21] permit only a single node to be 
released at a time by reducing the reaction force at the crack tip node 
over several steps in some prescribed manners when a fracture criterion, 
such as a critical crack tip opening displacement (CTOD) or crack tip 
opening angle (CTOA) , is reached. These methods require continual external 
monitoring and interruptions of the computations which become very time- 
consuming and inconvenient for simulation of crack growth over distances 
of ten or more elements. Hoff et al . [22] introduced a technique employ- 
ing spring and gap elements to circumvent this problem, the method 
requires a subroutine to control the openings of the gap elements and thus 
does not put crack propagation under total control of the external load- 
ings. In this study, we use nonlinear springs that have no restraining 
forces beyond a certain critical strain to imitate the failure charac- 
teristics of the cohesive material in the interlayer. Thus when a non- 
linear spring joining a node along the crack path to the foundation 
’fails', a node is released and the crack advances. In this way, crack 
propagations over distances of tens of elements (or more) are easily simu- 
lated requiring no external monitoring whatsoever. 
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The present problem is considered as a precursor to the more general 
one wherein a yield criterion based on a strain softening material is 
used. In that more general case one needs to be concerned with the highly 
triaxial stress state at the crack tip as well as with the consequences of 
locally unstable material behavior on global stability. In cases where 
the process zone is 'narrow' compared to the size of the specimen con- 
sidered, a 'boundary layer' type model in which the damage-softened 
material behavior is confined to a thin layer adjacent to the crack plane 
may prove to be a realistic and yet computationally inexpensive scheme for 
studying crack growth. Certain rubber-toughened materials seem to obey 
this model very closely. Much work remains before the appropriate contin- 
uum models for enigneering materials of interest are experimentally iden- 
tified and numerically verified under a wide range of loading conditions. 

In the next section, we discuss the problem of the beam on a non- 
linear foundation, its numerical solution schemes and the results for both 
stationary and propagating 2 3 cracks. The finite element model is presented 
in Section 3 and is followed by comparisons with the beam equation model 
in Section 4. Applications of the findings are demonstrated in Section 5, 
with conclusions summarized in Section 6. 


2. THE BEAM ON A NONLINEAR FOUNDATION MODEL 


The model and relevant geometric and material definitions are dep- 
icted in Figure 1. We start with the equation for a Bernoulli-Euler beam 
(with constant El) resting on a nonlinear foundation 


El 



q(w) 


0 


( 2 . 1 ) 


where E is the Young's modulus of the beam, I = ^bh 3 (for rectangular 

3. In the present investigation, no inertia (dynamic) effects are considered, 
therefore 'propagation' refers to 'quasi-static' propagation throughout. 
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cross section); b and h are the width and height of the beam, respec- 
tively, w is the vertical displacement of the beam neutral axis, and q(w) 
the nonlinear restoring force per unit length of the foundation. We con- 
sider q(w) such that, typically, the small strain response corresponds to 
an "elasticity" modulus that is approximately an order of magnitude 
smaller than that of the beam. 

Equation (2.1) is to be solved subject to boundary conditions at the 
loading end and at a distance far away from the loading end where the con- 
ditions are those of a semi-infinite beam on a linearly elastic founda- 
tion. This may always be assumed if the 'uncracked portion' of the beam 
is longer than two to three times the exponential decay length (l/j8; see 
equation (2.2)). 

2.1. Solution Schemes. 

The equation to be solved is nonlinear and involves boundary condi- 
tions at two points. A standard technique is the shooting method [23], 
which seeks the proper boundary conditions at one end point that also 
satisfy the boundary conditions at the other end point by Newton's itera- 
tion based on the Jacobian matrix formed by the miss-hits at the other end 
point and the current boundary conditions chosen at the first end point. 
The usual difficulty encountered is in obtaining convergence without 
requiring an excessive number of iterations; the rate of convergence 
depends strongly on the proximity of the initial guess of the boundary 
conditions at the first end point to the correct boundary conditions. 

In the present problem, we take advantage of the existence of the 
solution for the beam on a linear foundation by starting the integration 
from the region where the foundation response is linear and an analytical 
solution is available. We integrate numerically up to the crack tip 
beyond which the beam is free of surface tractions except for the very 
end. For this unloaded section a simple beam solution can again be 
exploited. For convenience we set up our coordinate system as shown in 
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Figure 2. 

The solution for a semi-infinite beam on an elastic foundation (i.e., 
q(w)=kw in (2.1)) with boundary conditions 

W(°°) = W ’ ( 00 ) = w' ' (00) = w ' ' ' (°°) = 0 


is of the form 


w(x) = e ^ x (Acos0x + Bsin0x) (2.2) 

1 

where 0 = 4 and A, B = constants. 

To start the integration of equation (2.1) from the elastic founda- 
tion region, we observe that the choice of the starting point is arbi- 
trary; for convenience we may, therefore, choose x to be zero such that 

w(0) = w q (2.3a) 

w'(0)= 0 (2.3b) 

which yields A = B = w q and 

w"(0) = -20 2 w q (2.3c) 

w' ' ' (0) = 40 3 w q (2.3d) 

Note that it is fastest to integrate from w^ < 0 into the nonlinear region 
(see Figure 2), and that w q has to be such that no yielding occurs in 
compression at the starting point x = 0. In this work, we assume that the 
yield strain in compression is equal to the yield strain in tension. It 
should be pointed out that by starting the integration from conditions 
derived from (2.2), we are guaranteed an exponentially decaying solution 
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as x -*■ + <»; this would not necessarily be the case if one were to start 
the integration from the loading end as any small error in the initial 
guess could cause the components of solution with terms proportional to 
e^ x to enter and make the solution unbounded as x * + ». Thus a poor ini- 
tial guess at the loading end could result in extremely slow convergence 
or no convergence at all. 

2.1.1. Stationary (Non-propagating) Case. For this case we seek the solu- 
tion for a geometry with a pre-cut crack of length 1 under increasing load 
P to the point where crack propagation is imminent. 


We obtain the solution by integrating 4 until § = -^ = 1 (with less 

r 5> 

than 0.1% error), here S is the shear for the free part of the beam, P is 
the end load, and M is the moment. We then obtain the rest of the solu- 
tion from (see Figure 2) 


3EI 


|W. |l 


+ w + 


(2.4) 


where w^ i S the displacement at the "crack tip"; w^. < w c , and | w ' t j is 
the slope of the neutral axis at the crack tip. When w = w , we reach 
the instant at which the crack is about to grow. 

2.1.2. Quasi-static Propagation Case. Here the scheme is to integrate 

along x until w = w £ (with less than 0.001% error) and then calculate the 
corresponding crack length and end conditions from: (see Figure 2) 


6 = 3EI + f w ’ c 1 1 + w c (2 ’ 5) 


4. The subroutine MODDEQ of the Caltech Computing Center which employs the 
Runge-Kutta-Gill method is used to solve a system of four first-order dif- 
ferential equations equivalent to equation (2.1). 
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where S c and M c are the shear and moment in the beam at the crack tip, and 

j w ' | is the slope of the beam there, 
c 

Other quantities such as the bending energy in the beam, , the 

total work done, W, the work done on the interlayer, W , the size of the 

c 

yield zone, a = |x c - x y | , and the 'secant' compliance, C s 6/p , are 
readily computed. Note that for a nonlinear system, such as the one we 

j * 

are dealing with, it is appropriate to define the compliance C as C = ; 

however in the present study, we intend to compare our results with previ- 
ous findings in other studies where C has always been defined as the 
secant compliance which, strictly speaking, is only correct for linear 
systems. As there is no obvious advantage in employing the more general 
definition here, the usual definition is used. 

2.2. Nondimensionalization. 


For presentation of the results it is useful to non-dimensional ize 
pertinent parameters. We choose the thickness ' d ' of the cohesive inter- 
layer as the natural length scale, i.e., define 


* 


w 


_ w 
“ d’ 



Equation (2.1) is then reduced to 5 


*( 4 ) 

w + CjW 


0 , 


(2.1a) 


with 


c ! • >*!ri is ! 3 


where E is the elastic (small strain) modulus of the cohesive foundation. 

L 


5. See reference 1 for detail. 
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We summarize here other relevant dimensionless quantities for future 
reference: 


6 * s 5/d, 1* s l/d, h* s h/d, 


b* 3 b/d 


* 

q 





* 

c 


* 


3 


* 

p 


and 7 * is the non-dimensionalized fracture energy (= area under the 

q-w curve) . 


2.3. Results and Discussions. 


We present next results for the stationary crack and the propagating 
crack. The material models considered for the interlayer are all 
piecewise-linear and possess the same fracture energy. The following set 
of data is used throughout: E is chosen as the Young's modulus for Alumi- 
num 2024 (10.6 x io 6 psi). E the small strain modulus of the cohesive 

C c 4c 

interlayer, is taken to be 1.41 * 10 psi. Also h = h/d = 10.0. These 
yield, = i,6 x io -4 . The material models are 7 shown in Figure 3. They 
are intended to examine the effects of various nonlinear material charac- 
teristics on the measurable macroscopic quantities. 

2.3.1. The Stationary Crack. In the following, the crack length is ten 

times the height of the beam. These results are presented in Figures 4-7. 

In Figure 4, the dashed lines exhibit the effects of the high compliance 

due to an initial low foundation resistance which subsequently increases 

* 

and becomes fairly constant as the yield zone, a , grows. (Thus the com- 
pliance decreases and then levels off.) These dashed portions are 

6. Note that E i s about two order of magnitude smaller than E. If E is 
greater than c a few percent of E, the beam may deform plastically near^the 
crack tip rendering the present beam equation model inadequate. 

7. The solution technique is capable of handling generally nonlinear q(w) , 
even though we only consider here piecewise-linear material models. 
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approximately extrapolated based on the compliance values from Figures 5. 
Clearly, the existence and the gradient of an softening tail (as in 
materials 2 & 3) can be easily identified in any of these five plots as 
the onset of crack growth is approached, i.e., as w w . This fact may 

L 0 

be used to approximately characterize q(w) from experimental data. In 

general, the compliance, C*, the slope at the crack tip, -w*' , and the 

* 

size of the yield zone, <x , are all monotonically increasing functions of 

* 


We shall discuss the application of this finding in connection with 
the characterization of q(w) in more detail in Section 5.2. 

2.3.2. The Quasi-statically Propagating Crack. Before dealing with the 
nonlinear case we consider first certain results for the idealized situa- 
tion where neither the interlayer nor the elasticity of the beam allows 
rotation of the 'built-in' end. We first derive the equations mentioned 
in Section 1 for the built-in beam case. Here the coordinate system is 
such that the load P is applied at x = 0, w(o) = 5 and w(l) = w'(l) = 0. 
We have for this "idealized" case, 


j = El! 

3EI 


and for the bending energy in the beam U , 


IT = } Mi^dx = -Ei- } x 2 dx 
b { 2EI ax 2EI 1 


> 2 1 3 


6EI 


( 2 . 6 ) 


(2.7) 


Griffith Condition : The energy release rate, G, is obtained from the 
potential energy n of the system (both beams of the DCB specimen con- 
sidered and hence the factor two) as follows 


G 


an = __a_ 
ai ai 


2[U b -P3] 



• P 2 ! 3 P 2 ! 2 

6EI ' El 


27 b 


Hence the moment at the crack tip required for crack growth is 

L 
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M c = pi = \|2El7b 

Combining (2.6) and (2.8), one finds 

3P 2 « = (2Yb) 3/2 \|ET. 
Rewriting (2.9), one obtains 


27 


1 I 13E_* ) 

b { El 


2 1 
- > 
J 


1/3 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


In order to study how a certain quantity changes as the crack pro- 
pagates quasi-statically, we denote the derivative with respect to crack 

length 1 by -j f 1 q3 . Note that P is not constant as 1 increases but 
decreases such that quasi-static crack growth is maintained. Adopting 
this convention and using (2.7) and (2.8), there follows for the bending 
energy of the beam (here only half of the DCB specimen is considered) 

M^l 

U b = 6EI = { V )lf 


whence 


dl 



( 2 . 11 ) 


For this ideal case, the rate of work done against interface forces as the 
crack propagates is 


dW c 

dl ! qs = ~ b (2-i2) 


where denotes the energy dissipated in the interlayer. Let W be the 

5 

total work done by P, i.e., W = / Pd6 , or dW = Pdi. Then in view of 

o 


(2.6) and (2.8) 
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d5, 

dl'qs 


!V 

3EI 


so that 


a, .pdt, .3£i 

dl'qs ‘dl'qs 3EI 



In view of (2.8), 


dW| m *y b 
dl'qs 3 7 


Combining (2.11), (2.12), and (2.13) yields 


dU h 

— | 
di qs 

dw c 

df'qs 


1 dW| 

4 dl'qs 

3 dW, 

4 dl 'qs 


(2.13) 


(2.14) 


Thus for the idealized (built-in) case, 25% of the total work done is 
stored in the beam as bending energy while 75* is dissipated as fracture 
energy at the interface. Note in passing that in this idealized case the 
stress intensity factor is 


K = (2\|3h 3/2 )Pl = (2\|3h 3/2 )M, (2.15) 

But by (2.8), M c = f( 7 ,EI), therefore Kj at fracture is a constant as is G 
(the energy release rate). 


Having considered this ideal case we consider now the situation where 
the nonlinear foundation allows both the displacement and the slope at the 
crack tip to be different from zero. The numerical results 8 in Figure 8 
show that the P-6 relation depends only on y and El as suggested by 

8. Here material #2 was not investigated, as we did not expect the use of its 
characteristic to contribute any additional insight into the problem. 
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$ 

equation (2.9) derived for the idealized case. Figure 9 indicates that M 

0 

* 

depends weakly on w and on the existence of a strain softening tail. The 

* i 

slope at the crack tip, -w , increases and rapidly approaches an asymp- 

* 

totic value as the crack grows. The yield zone size a is shown to 

* * Q * *1, 

decrease and becomes constant as 1 /h increases. In Figure 10, P 5 74 is 

* * 

shown to approach a constant value as 1 /h increases. In view of equation 
(2.9) and the behavior of P*6*^ exhibited in Figure 10, it seems possible 
that the fracture energy y in the case where the beam supports rotates may 
still be approximated by equations of the same form as (2.10) and (2.10a). 
With this in mind, define 


27 


= i { cagll).!} 

- b 1 EI j 


1/3 


(2.16) 



s JL ( 3 P*V ) 2 / 3 


(2.16a) 


~ * 

In Figure 14, we compare the value of y , the fracture (surface) energy as 

* 

computed from equation (2.16a), to the exact value y as calculated from 
the area under the q-w curve of the input data. Note that the vertical 
ordinates in Figures 13 and 14 encompass very narrow ranges of values. 10 


In these non-ideal cases, it is found that equations (2. 9) -(2. 14) 
remain valid with small errors which decrease rapidly with increasing 

crack length. The reason for this behavior is that the crack propagation 

* 

has reached an asymptotic condition in the sense that neither a nor -w is 

c 

changing much. To see this consider the following argument: 


For the idealized case, equation (2.10) can be written in functional 
form as 

9. For quantitative display, see reference 1. 

10. The accuracy of these results is demonstrated quantitatively in reference 

22 . 
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F(P, 5, 7 , El) = 0 (2.17) 

Note that, the crack length 1 does not appear explicitly as it is already 
fixed through the equilibrium condition by specifying both P and 6 . In 
the event that the support rotates (the nonlinear foundation case), the 
only additional physical parameter 11 which could be involved in (2.17) is 

' dw c, 

w c> however << l for 1/ h > 5 (see Figure 11). Thus, the slope at 

the crack tip hardly changes as the crack propagates and, therefore, 

(2.17) ' still holds with good approximation for 1/ h > 5. 

One could also think of the built-in case as a special case where 

I 

q(w) = 76 (w) = a delta function of integral measure 7 and w = a = 0 . 

c 

Thus, for general q(w), one only needs to reach the asymptotic conditions 

t 

for a long crack when w and a are approximately constant for equation 

(2.17) to hold. 


3. FINITE ELEMENT MODEL 


The finite element program ABAQUS (version 4) running on a VAX-11/780 
was used to model the beam and the nonlinear foundation. We use 4-noded 
bilinear plane elements to model the beam, and nonlinear springs restrain- 
ing vertical motion of the beam to model the foundation. The 4-noded ele- 
ment is chosen for convenience in discretizing the foundation and in 
interpolating to locate the crack tip. The force-displacement relations 
of the nonlinear springs are patterned precisely after the material models 
shown in Figure 3. 


11. Note that we may disregard w as an additional parameter since it is a 
constant for a given q(w) under consideration. It should also be pointed 
out that in all cases considered |w^| is very small (< 0.014) i.e. less 
than 0.80° in rotation at the crack tip, but the contribution of this 
small slope to the end displacement 5 is not at all negligible, for 
instance, by equation (2.5) for 1/ h = 10, h/ d = 10; |wjl = 0.36. 
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The discretization of the foundation is accomplished by attaching two 
identical nonlinear springs with total restraining forces equal to that of 
a continuous foundation to the bottom two nodes of each element of the 
beam at the interface. (For 8-noded biquadratic elements, this simple 
scheme would not work.) 

The three meshes used for the beam are shown in Figure 12. Mesh #1 
has 402 elements, 67 springs, 476 nodes, and a total of 952 degrees of 
freedom with the smallest element size of O.lh * 0.05h. Mesh #2 has 540 
elements, 90 springs, 637 nodes, and 1274 degrees of freedom with the same 
smallest element size as Mesh #1. Mesh #3 has 480 elements, 120 springs, 
605 nodes and 1210 degrees of freedom with smallest element size of 0.25h 
x O.lh. (For Mesh #2, an iteration takes approximately 75 CPU seconds and 
each increment requires 2 or 3 iterations.) 

Mesh #1 is designed to capture the displacement details around the 
* * 

crack tip for i / h =10.0 while mesh #2 is intended to study a shorter 

4c 4c 

specimen and shorter crack lengths, i.e., 1 / h <3.5. Mesh #3 is a more 

uniform mesh than the previous two and gives better results for a larger 
* * 

range of 1 / h (from 0.0 to 10.5). 

The loading is achieved by prescribing the end displacement (5) at 
the first node on the bottom line of the beam; the reaction force is then 
equal to P. The crack length is obtained by linear interpolation between 
a node where w > w^’ and an adjacent node where w < w c . (Tills convenience 
would be lost if 8-noded biquadratic elements were used.) Convergence is 
considered attained when the forces at all nodes except those with 
prescribed displacements fall below 0.2% of the typical applied force 
values (in this case P). Since the beam is linearly elastic and the non- 
linear springs simply supply the proper boundary conditions, the typical 
residuals for the convergent solutions are = 10 p which were much 
smaller than the set tolerance. 
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4. COMPARISONS OF THE BEAM MODEL WITH FINITE ELEMENT RESULTS 

In this section, we examine the effects of including shear deforma- 
tions through the finite element model while neglecting them in the beam 
equation model. The finite element results show that, excepting the beam 
end where the point force is applied, the stress components and u 22 in 
the beam are always less than 10% of the maximum value of <r The max- 

imum equivalent (von Mises) stress in the beam only reaches 53% of the 
yield stress (for Aluminum 2024, <J y = 50. * 10 3 psi). For beam materi- 
als with a low yield stress one must ensure that no yielding occurs in 
order for the present analysis to be applicable. 

All comparisons shown here are for material model 3 (see Figure 3). 
First, the displacement contours of the top, center and bottom of the beam 
are compared with the solution to the beam equation (the neutral axis or 
center line displacement) at the same crack length to determine whether 
measurments can be taken at the top of a DCB specimen instead of at the 
bottom (interlayer) where it would be more difficult to measure; Figures 
13 gives the percent differences of the two results. It is seen that the 
calculated displacements at the top, center or bottom of the beam are vir- 
tually identical; this is true as long as plastic deformations are absent. 
Thus the beam equation model yields results within a few percent of the 
finite element analysis for most of the length of the beam except at loca- 

* * 5 ^ 

tions ahead of the crack tip. In that region (x / h > 11.0) w is very 
small and the mesh for the finite element analysis becomes coarse thus 
causing the percent error to rise. 

* * 

The P -6 relationships obtained from both methods are compared in 

Figure 14. The finite element results for Mesh #2 and #3 give also the 
▼ * 

P -6 points before crack propagation commences ( the peak in the upper 
left hand corner) . The best agreement is obtained by using a uniform mesh 
(#3). 
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5. APPLICATIONS 


We discuss next two practical applications of the findings, one is a 
rational method for determining the fracture surface energy, the other 
which is the primary goal of this study, namely, the rough characteriza- 
tion of q(w) . 

5.1. Determination of the Surface Energy. 

In order to assess the delamination strength of composites or the 
strength of adhesives, it is desirable to determine the energy expended in 
causing crack propagation per unit of failure surface. Because of its 
relative simplicity the cleavage or double-cantilever-beam (DCB) specimen 
is often employed for this purpose. It has long been recognized that the 
simple formula derived on the basis of a built-in condition at the crack 
tip (see Section 2.3.2 and equation (2.16)) is not necessarily applicable 
when the notation (slope) at the crack tip is nonzero. To circumvent this 
difficulty, Berry [2] introduced an ad hoc power-law scheme which is pur- 
ported to be more accurate than the equation (2.16). However, as the 
results for the propagating crack in Section 2.3.2 suggest (see Figure 
14), equation (2.16) can be used to calculate y (= y ) conveniently and 
accurately. 

We believe that the application of equation (2.16) is a more rational 
method for the determination of the surface energy than Berry's scheme 
which has been widely accepted. In the following, we will show numeri- 
cally that Berry's scheme is an approximation to the results obtained 
through the solution of the beam equation with the finite rotation at the 
crack tip incorporated through the specification of q(w). 12 We start by 
outlining both schemes. 

12. Note that regardless of the detailed characteristics of q(w), a finite 
rotation (slope) always exists at the crack tip and that this rotation 
contributes substantially to the end displacement as indicated in footnote 
11 . 
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Metho d (A): Berry [2] asssumes that the compliance is a power-law function 
of the crack length, 1 , in the form 


C 



N 


where a depends on El. Upon plotting log C versus log 1 and determining 
the slope N, one finds the energy release rate G to be 



Let us adopt the notations (superscripts refer to method A) 

G (A) = 2 y (A) b = 


Thus one needs to measure several sets of P,5 and 1 to be able to 
determine N resonably accurately from the plot of log C versus log 1 . 


Method (B) : The method proposed here is the application of equation 

(2. 16), 13 i.e., we let, 


g (B) 


27 


(B) v 



(5.2) 


where 7 ^ = 7 defined by (2.16). Both schemes are applied to experimen- 
tal data on composite delamination taken from reference 26 and are com- 
pared in Table I . 


13. This method was proposed and used earlier in [24,25] for cases where the 
'built-in' conditions are assumed to be valid. 
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Table I 

7 

Material : Unidirectional Composite T300/5208, E = 1.9 * 10 psi. 


N = 2.68, 

(N is 

obtained 

from 

the plot of log C 

versus 

log 1). 


h = 0.119 

in. , 

b = 0.205 

in. , 

El = 547.0 lb-in 2 




1/h 

P(lb) 

l(in) 


5 ( in) 

2*y( A ) 

27 (B) 

2v tA) 

2 7 < B > 

22.3 

3.05 

2.65 


0.105 

1.58 

1.22 

1.30 

24.9 

3.12 

2.96 


0.135 

1.86 

1.49 

1.25 

33.8 

2.31 

4.02 


0.220 

1.65 

1.38 

1.20 

44.9 

1.75 

5.34 


0.375 

1.61 

1.36 

1.18 

47.6 

1.82 

5.67 


0.425 

1.78 

1.56 

1.14 

52.5 

1.50 

6.25 


0.480 

1.51 

1.31 

1.15 

56.9 

1.42 

6.77 


0.540 

1.48 

1.31 

1 . 13 

62.8 

1.37 

7.47 


0.690 

1.65 

1.47 

1.12 





2V ave (lb/in) 
Standard 

1.65 

1.39 






Deviation 

0.12 

0.10 



One notes from Table I that the results for y based on 3erry's scheme 
(Method A) are consistently higher than for the scheme proposed here 
(Method B) by roughly 20% on the average. The question now arises whether 
one method yields a better estimate for y than the other. 

Without knowing from another independent .source a very reliable value 

of y, we resort to an indirect argument based on the numerical solutions 

for the model of a beam on a nonlinear foundation. To do this, we calcu- 

* * 

late N numerically as a function of 1 /h by following Method (A). -- If 

Berry's assumption is correct, N should be nearly constant. As depicted 

* * 

in Figure 15, N increases with 1 /h and approaches N = 3 which 


14. The values of N are the local slopes of the log C versus log 1 plot. 
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4c * 

corresponds to the idealized case (see equation (2.6)) as 1 /h -► «>. it 
is evident that N varies considerably for all three material models stu- 
died especially for shorter crack lengths. Therefore the value of N 
obtained using Berry's scheme is simply an approximation. 

The fact that we find N in Berry's method to be a function of 1 indi- 
cates that the power law assumption is not really acceptable (since an 
analytical basis is lacking in the first place) and that the proposed 
method (B), allowing for the rotation of the 'beam end', is inherently 
more rational.*® 

r (B) 

To further confirm the above point, we list values of u = in 

G ( A ) 27 ( A ) 

the last column of Table I. The surprisingly monotonic decreasing trend 
* 

with increasing 1 of the above ratio in the data suggests that a high 
value for N (N = 2.68) was obtained from the log C versus log 1 plot. 
This high value for N corresponds to a rather long crack length (see Fig- 
ure 15), and would thus yield an over estimate for 7 for shorter crack 
length 1 since N is assumed to be constant in equation (5.1). This above 
discussion explains the decreasing trend in the ratio in the last column, 
since for longer crack lengths 1, this high value N becomes closer to the 
proper local values of N's that would theoretically yield the correct 
value for 7.*® 


15. It turns out that even if one uses a continuously varying N for calculat- 
ing 7 through equation (5.1), the results are still less accurate than the 
simple application of equation (5.2). See reference 1 for quantitative 
details . 

16. This comparison was made on two sets of experiments which yielded the same 
qualitative results, see reference 1 for more detail. The experiments 
from which the data are taken were carried out on nonsymmetric DCB speci- 
ments, i.e., the heights of the two beams were not always equal; we assume 
5 to be the average of the sum of the 5's of each beam. Details are given 
in reference 26. 
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5.2. Characterization of q(w). 

We next discuss how the nonlinear characteristics of q(w) may be 

estimated based on information gathered from basic material testings and 

tests on DCB specimens. The properties of q(w) which are of interest are 

E w , w , 7 as well as the existence and gradient of an unloading tail, 
v y l 

We outline below the methodology by which the characterization may 
proceed: 

a. E ^ and w can be determined from uniaxial tests. 

c y 

b. w c may be obtained from ultimate strain tension tests or optical 
measurements of the displacement contour from the DCB specimen 
tests . 

c. 7 is easily computed using equation (2.16) or equation (5.2) as dis- 
cussed at length in Section 5.1. 

d. The existence and the gradient of an unloading tail can be identified 

and estimated by plotting C, -w t and a as a function of P using 

DCB specimen test results for the stationary crack case (see Section 
2.3.1, Figures 5-7 ) . 

These determinations are sufficient to provide a bound on the shape 

of the function q(w) . Refinements of this approximate characterization 

1 ft 

can then be obtained by solving the beam equation for various assumed 
q(w)’s (within the above bound) to better match experimental results. 


17. Note that from results depicted in Figure 16, the displacement profile can 
be conveniently taken at the top surface of the beam since it is virtually 
identical to the profiles at the center and bottom surface of the beam (as 
long as there is no plastic deformation). 

18. We mention in passing that the solutions of the beam equation are rela- 
tively inexpensive compared to the finite element solutions. 
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6. CONCLUDING REMARKS 

We have established a possible experimental-analytical procedure by 
which certain nonlinear (cohesive) material characteristics may be deter- 
mined approximately. In addition, a rational method for determining the 
surface (or fracture) energy in DCB specimen tests which accounts for the 
rotation at the crack tip ('built-in' end) and supplants Berry's method is 
proposed . 
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9. FIGURE 

Figure 1. 
Figure 2. 
Figure 3. 
Figure 4. 

Figure 5. 
Figure 6. 
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CAPTIONS 


Problem definition of the geometry and material parameters. 


Solution schemes. 


Material models (1, 2, 3 and 4). 


End load versus end displacement for the stationary crack 
(material models 1, 2, 3, and 4 as indicated). For explana- 
tion of the dashed portions of the P-6 curves see text. 


Compliance as a function of the end load for the stationary 
crack (material models 1, 2, 3 and 4 as indicated). 


Crack tip displacement as a function of the end load for the 
stationary crack (material models 1, 2, 3 and 4 as indicated). 


Slope at the crack tip as a function of the end load for the 
stationary crack (material models 1, 2, 3 and 4 as indicated). 


End load versus end displacement for the propagating crack 
(material models 1, 3 and 4 as indicated). Note that the 

relation displayed is not affected by the detailed charac- 
teristics of q(w) . 


Moment at the crack tip as a function of crack length for the 
propagating crack (material models 1, 3 and 4 as indicated). 


3jC 

P 6 as a function of crack length for the propagating crack 
(material models 1, 3 and 4 as indicated). 


~ * * 

7/7 as a function of crack length for the propagating crack 
(Material models 1, 3 and 4 as indicated) 


Figure 12. Finite element meshes. 
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Figure 13. Comparison between Finite Element and Beam Equation: Percent 
difference in vertical displacement as a function of x (for 
the same crack length). 


Figure 14. Comparison between Finite Element and Beam Equation: End load 
versus end displacement for 3*different Finite Element meshes. 
(The higher discrepancy at 5 /h >2.7 is due to the coarse 
mesh region.) 


Figure 15. The exponent N in the assumed power law C = al N as a function 
of crack length for the propagating crack (material models 1, 
3 and 4 as indicated). 



PROBLEM DEFINITIONS 



E = Beam Modulus 

E c = Cohesive Elastic Modulus 

p = End Load 

5 = End Displacement 

w = Vertical Displacement 

6 - Crack Length 
L = Beam Length 

a = Yield Zone Size 


h = Beam Height (thickness) 
b = Beam Width 
d = Cohesive Layer Thickness 
q(w)= Cohesive Restoring Force 
Wy = = Yield Strain 

w*y = ■^ L = Unloading Strain 

# w 

w c = - Critical (crack) Strain 


Figure 1. Problem definition of the geometry and material parameters. 
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Figure 2. Solution schemes. 


(All models have the same y 




Figure 3. Material models (1, 




Figure 4. End load versus end displacement for the stationary crack 
(material models 1, 2, 3, and 4 as indicated). For explana- 
tion of the dashed portions of the P-8 curves see text. 







Figure 6. Crack tip displacement as a function of the end load for the 
stationary crack (material models 1, 2, 3 and 4 as indicated). 




Figure 7. Slope at the crack tip as a function of the end load for the 
stationary crack (material models 1, 2, 3 and 4 as indicated). 






Figure 8. End load versus end displacement for the propagating crack 
(material models 1, 3 and 4 as indicated). Note that the 

relation displayed is not affected by the detailed charac- 






Figure 10. P 6 as a function of crack length for the propagating crack 
(material models 1, 3 and 4 as indicated). 




Figure 11. 7/7 as a function of crack length for the propagating crack 

(Material models 1, 3 and 4 as indicated) 



Beam on Nonlinear Cohesive Foundation 

Finite Element Meshes 





Figure 12. 


Finite element meshes. 




Figure 13. Comparison between Finite Element and Beam Equation: Percent 
difference in vertical displacement as a function of x (for 
the same crack length). 
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Figure 14. Comparison between Finite Element and Beam Equation: End load 
versus end displacement for 3 # di£ferent Finite Element meshes. 
(The higher discrepancy at 8 /h >2.7 is due to the coarse 
mesh region.) 
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Figure 15. The exponent N in the assumed power law C = al N as a function 
of crack length for the propagating crack (material models 1, 
3 and 4 as indicated). 



